[velP,velM,velP_ind,velM_ind,velP_raw,velM_raw]=get_vel_ind_from_adata(aux_data);


%% treebagger stuff
cnt=0;
for knd=1:4
    for ind=1:size(ROIs{knd},2)
        cnt=cnt+1;
        aaa(:,cnt)=calc_dFF(ROIs{knd}(ind).activity);
    end
end


%aaass=aaa(:,round(rand(ind,1)*93+1));
fprintf('\n000')
for ind=1:size(aaa,2)
    fprintf('\b\b\b\b%4u',ind)
    tp=aaa(:,ind);
    mtp=median(tp);
    stdtp=std(tp);
    tp=psmooth(tp);
    tp=tp/mtp;
    %tp(tp<mtp+stdtp)=mtp;
    aaass(:,ind)=tp;
end

bbb=[];
bbb(:,1)=smooth2(velM(frame_times(2:4:end-4)),10);
% bbb(:,2)=smooth2(velM(frame_times(2:4:end-4)),10);

B=TreeBagger(32,aaa(10000:15000,:),bbb(10000:15000),'method','regression');
ccc=predict(B,aaa);

[r,p]=corrcoef(bbb(1:5000),ccc(1:5000));
fract_var_expl2(ind,knd)=r(2)^2;

figure;
plot(ccc)
hold on
plot(bbb,'r')


bbb=[];
bbb(:,1)=smooth2(velM(frame_times(2:4:end-4)),10);
%bbb(:,2)=smooth2(velM(frame_times(2:4:end-4)),10);
B=TreeBagger(32,bbb(5001:10000),aaa(5001:10000,20),'method','regression');
ccc=predict(B,bbb);

figure;
plot(aaa(:,20))
hold on
plot(ccc,'r')


bbb=[];
bbb(:,1)=smooth2(velM(frame_times(2:4:end-4)),10);
%bbb(:,2)=smooth2(velM(frame_times(2:4:end-4)),10);
B=TreeBagger(32,bbb(10001:15000),aaa(10001:15000,20),'method','regression');
ccc=predict(B,bbb);

figure;
plot(aaa(:,20))
hold on
plot(ccc,'r')

%% population response to fb mismatch
bbb2=[];
bbb2(:,1)=smooth2(velM(frame_times(2:4:end-4)),10);

pertub = (aux_data(3,frame_times(2:4:end-4))')>0.5;  
%%% taking layer 2 as reference is enough for now

pertub_on = strfind(pertub(1:5000)',[0,1]);

pertub_response = zeros(41,size(aaa,2)*length(pertub_on));
for kknd = 1:length(pertub_on)
    pertub_response(:,1+size(aaa,2)*(kknd-1):size(aaa,2)*kknd) = aaa(pertub_on(kknd)-10:pertub_on(kknd)+30,:);
end

 figure,imagesc(pertub_response')
 figure,plot(mean(pertub_response,2))  
 temp = reshape(pertub_response,size(pertub_response,1),size(pertub_response,2)/length(pertub_on),length(pertub_on)); 
figure,imagesc(mean(temp,3)')
% figure,plot(mean(mean(temp,3),2))

sss=mean(pertub_response,2);
figure;plot(sss-mean(sss(1:10)))
ylim([-0.01 0.04])

%% corrected for pertubations only during running

aaa = [];
cnt=0;
for knd=1:4
    for ind=1:size(ROIs{knd},2)
        cnt=cnt+1;
        aaa(:,cnt)=psmooth(calc_dFF(ROIs{knd}(ind).activity));
    end
end
[velP,velM,velP_ind,velM_ind,velP_raw,velM_raw]=get_vel_ind_from_adata(aux_data);

bbb2=[];
bbb2(:,1)=smooth2(velM(frame_times(2:4:end-4)),10);

pertub = (aux_data(3,frame_times(2:4:end-4))')>0.5;  
%%% taking layer 2 as reference is enough for now

pertub_on = strfind(pertub(1:5000)',[0,1]);

pertub_response = zeros(41,size(aaa,2)*length(pertub_on));
cnt = 1;
for kknd = 1:length(pertub_on)
    if sum(bbb2(pertub_on(kknd)-10:pertub_on(kknd)+10) > 0.002) > 16   %% running threshold and no. frames 
        pertub_response(:,1+size(aaa,2)*(cnt-1):size(aaa,2)*cnt) = ...
            aaa(pertub_on(kknd)-10:pertub_on(kknd)+30,:);
        cnt = cnt + 1;
    end
end
cnt = cnt - 1
pertub_response = pertub_response(:,1:size(aaa,2)*cnt);

sss=mean(pertub_response,2);
figure;plot((sss-mean(sss(1:10)))*100)
ylim([-2 4])
set(gca,'XTickLabel',[-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('seconds')
ylabel('dF/F [%]')
title('population response to feedback mismatch onset')

 temp = reshape(pertub_response,size(pertub_response,1),size(pertub_response,2)/cnt,cnt); 
 figure,imagesc(mean(temp,3)')
 figure,plot(mean(temp,3))
 figure,plot(mean(mean(temp,3)'))
 pertub_rois_mean = mean(temp,3);
 
 %% running onsets

run_onset = bbb2(1:5000)>0.002;
a1 = [1,run_onset',1];
stloc = strfind(a1,[1 0]);
endloc = strfind(a1,[0 1]);
countseqs = endloc - stloc;
run_onset(stloc(countseqs == 1)) = 1;

a1 = [0,run_onset',0];
stloc = strfind(a1,[0 1]);
endloc = strfind(a1,[1 0]);
countseqs = endloc - stloc;
run_onset(stloc(countseqs == 1)) = 0;

run_onset = strfind(run_onset',[0 1]);
run_response = zeros(41,size(aaa,2)*length(run_onset));
cnt = 1;
for kknd = 1:length(run_onset)
%     if sum(bbb2(run_onset(kknd)-10:run_onset(kknd)+10) > 0.001) > 16   %% running threshold and no. frames 
        run_response(:,1+size(aaa,2)*(cnt-1):size(aaa,2)*cnt) = ...
            aaa(run_onset(kknd)-10:run_onset(kknd)+30,:);
        cnt = cnt + 1;
%     end
end
cnt = cnt - 1
run_response = run_response(:,1:size(aaa,2)*cnt);
sss2=mean(run_response,2);
figure;plot((sss2-mean(sss2(1:10)))*100)
hold on 
plot((sss-mean(sss(1:10)))*100,'g')
ylim([-2 4])
set(gca,'XTickLabel',[-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('seconds')
ylabel('dF/F [%]')
title('population response to feedback mismatch onset and running onset')

temp = reshape(run_response,size(run_response,1),size(run_response,2)/cnt,cnt); 
running_rois_mean = mean(temp,3);

 
 %% corr velP and velM to activity
 
ddd2 = [];
ddd2(:,1) = smooth2(velP(frame_times(2:4:end-4)),10);

pb_ind = 5001:10000;
corr_result = zeros(size(aaa,2),2);
fprintf('\n000')
for ind = 1:size(aaa,2)
    fprintf('\b\b\b\b%4u', ind)
    corr_result(ind, 1) = corr2(aaa(pb_ind, ind), bbb2(pb_ind));
    corr_result(ind, 2) = corr2(aaa(pb_ind, ind), ddd2(pb_ind));
end

figure,plot(corr_result(:,1),corr_result(:,2),'.')
ylim([-1 1]); xlim([-1 1])
axis square
xlabel('r_{running vs. activity}')
ylabel('r_{visual flow vs. activity}')

%% summmary
regions = {};
rois_per_region = [];
pertub_rois_mean_all = [];
running_rois_mean_all = [];
corr_result_all = [];

rois_per_region = [rois_per_region size(pertub_rois_mean,2)];
pertub_rois_mean_all = [pertub_rois_mean_all pertub_rois_mean];
running_rois_mean_all = [running_rois_mean_all running_rois_mean];
corr_result_all = [corr_result_all; corr_result];
regions{end+1} = fnames{1}; 

save('c:\temp\axons_ana_12115smooth.mat','rois_per_region','pertub_rois_mean_all','running_rois_mean_all','corr_result_all','regions');




figure,plot(corr_result_all(:,1),corr_result_all(:,2),'.')
ylim([-1 1]); xlim([-1 1])
axis square
xlabel('r_{running vs. activity}')
ylabel('r_{visual flow vs. activity}')
set(gca,'XTick',-1:0.5:1,'YTick',-1:0.5:1)


sss3=mean(mean(running_rois_mean_all,3)');
sss4=mean(mean(pertub_rois_mean_all,3)');
figure;plot((sss3-mean(sss3(1:10)))*100)
hold on 
plot((sss4-mean(sss4(1:10)))*100,'g')
ylim([-2 4])
set(gca,'XTickLabel',[-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('seconds')
ylabel('dF/F [%]')
title('population response to feedback mismatch onset and running onset')

%% plots
aaa_calc = [];
for ind = 1:length(ROIs{2})
    aaa_calc(ind,:) = calc_dFF(ROIs{2}(ind).activity);
end
figure,subplot(3,1,1:2);imagesc(aaa_calc)
set(gca,'clim',[0.5 3])
subplot(313);plot(smooth2([ddd2(1:10000)/max(ddd2); zeros(5000,1)]',100)-1.1,'r')
hold on
subplot(313);plot(smooth2((bbb2/max(bbb2))',100),'b')
subplot(313);plot(pertub_on,ones(size(pertub_on))-2.3,'k.')
xlim([1 15000])

cell_no = 64;
figure, plot(smooth2(psmooth(aaa_calc(cell_no,:))',10)-1,'g')
hold on
plot(smooth2((bbb2/max(bbb2))',100)-1.1,'b')
plot(smooth2([ddd2(1:10000)/max(ddd2); zeros(5000,1)]',100)-2.2,'r')
plot(pertub_on,ones(size(pertub_on))-3.3,'k.')
xlim([1 15000])

figure,imagesc(ROIs2image(ROIs{2}(cell_no),size(template{2}),0))
colormap gray
axis equal 
axis off

